#!/usr/bin/R
	
#begin log file
	path='/pondSim/'
	paramDir='/pondSim/params/'
	inputDir='/pondSim/input/'
	logfile='/pondSim/pondSim.log'
	write(file=logfile,x=paste('PONDSIM: ',date()))

#read model parameters###########################################################
	write(file=logfile,append=T,paste('READING INPUT/PARAMETERS',date()))
	source(paste(path,'pondSim.module.R',sep=''))
	read.params()
	read.pumps()
	read.input()

#set tolerances##################################################################
	tol.1=5e-3 #max change in ft during single timestep (dictates max delta.t)
	tol.2=1e-3 #min change in ft during single timestep (dictates condition where delta.t increased)
	tol.3=1 #dictates min delta.t
	t.factor=1.2 #factor by which delta.t increased/decreased
	
#intialize variables#############################################################
	write(file=logfile,append=TRUE,x=paste('ALLOCATING MEMORY',date()))
	n.times=5*24*60*60 #max number of timesteps (5 days)
	total.t=array(0) #total time simulated (sec)
	h=array(0,dim=c(n.ponds,n.times)) #water surface elev above MSL
	pond.inflow=pond.outflow=array(0,dim=c(n.ponds,n.times)) #net inflow to each pond (cfs)
	ff=array(0,dim=c(n.cons,n.times)) #flow through control structures
	K=array(0,dim=c(3,n.ponds)) #dummy variables for RK routine
	h[,1]=p.surf #set initial conditions
		
#perform simulation###############################################################
	write(file=logfile,append=TRUE,x=paste('BEGIN SIMULATION',date()))
	t=day.old=0
	while(TRUE) {	
		t=t+1
		day=1+as.integer(total.t[t]/24/60/60)
		hour=1+as.integer(total.t[t]/60/60)
		quarter.hour=1+as.integer(total.t[t]/60/15)
		h[p.bay,t]=SSE[(day-1)*24*4+quarter.hour] #sea surface elev

		#Check for end of simulation################################################
                if (total.t[t]>=t.horizon | hour==(24*5)) {
                        write(file=paste(inputDir,'pond.surf',sep=''),ncolumns=1,h[,t-1])
                        update.data()
                        write(file=logfile,append=T,x=paste('SIMULATION SUCCCESSFUL',date()))
                        print('SIMULATION SUCCESSFUL')
                        make.bargraph()
                        break
                }

		#make daily adjustment to v-notch based on target flowrate to allen ('ff.allen')
		if (day>day.old) {
			if (day==1) con.crest[con.allen]=cl.flow.to.crest(plant.inflow[day])
			else con.crest[con.allen]=cl.flow.to.crest((plant.inflow[day]+plant.inflow.old)/2)
			plant.inflow.old=plant.inflow[day]
			day.old=day
		}
	
		#RK STEP 1
		foo=find.flow(h[,t])
			pond.inflow=foo[[1]]+(PP[day]*p.area)
			pond.outflow=foo[[2]]+(ET[day]*p.area)
			pond.inflow[1]=pond.inflow[1]+plant.inflow[day] #add plant inflow to first pond
			K[1,]=delta.t/p.area*(pond.inflow-pond.outflow)
			h.temp=h[,t]+(1/3)*K[1,]
		#RK STEP 2
		foo=find.flow(h.temp)
                        pond.inflow=foo[[1]]+(PP[day]*p.area)
                        pond.outflow=foo[[2]]+(ET[day]*p.area)
                        pond.inflow[1]=pond.inflow[1]+plant.inflow[day] #add plant inflow to first pond
                        K[2,]=delta.t/p.area*(pond.inflow-pond.outflow)
                        h.temp=h[,t]+(2/3)*K[2,]
		#RK STEP 3
		foo=find.flow(h.temp)
                        pond.inflow=foo[[1]]+(PP[day]*p.area)
                        pond.outflow=foo[[2]]+(ET[day]*p.area)
                        pond.inflow[1]=pond.inflow[1]+plant.inflow[day] #add plant inflow to first pond
                        K[3,]=delta.t/p.area*(pond.inflow-pond.outflow)
                        h.temp=h[,t]+(1/4)*K[1,]+(3/4)*K[3,]
		
		#Implement variable timestep###############################################
		#if heads change by alot, reduce delta.t
		if (any(abs(h.temp-h[,t])>=tol.1) & (delta.t/t.factor)>tol.3) {
			if ((delta.t/t.factor)>tol.3) {
				delta.t=delta.t/t.factor 
				t=t-1
				next
			}
			else {
				write(file=logfile,append=TRUE,
					x=paste('min timestep not small enough to capture change in levels',date()))			
				write(file=logfile,append=TRUE,x=paste('SIMULATION STOPPED',date()))
				stop_()
			}
		}
		#if heads change by a small amount, increase delta.t
		else if (all(abs(h.temp-h[,t])<=tol.2)) {
			delta.t=delta.t*t.factor
		}
		
		#if timestep ok, store new head value#######################################
                h[,t+1]=h.temp
                total.t[t+1]=total.t[t]+delta.t

		#Check for high/low water levels############################################
		if (any(h.temp<p.min.surf)) { 
			write(file=logfile,append=T,x=paste('low head in ponds:',p.names[which(h.temp<p.min.surf)]))
			write(file=logfile,append=T,x=paste('min head=',p.min.surf[which(h.temp<p.min.surf)],
				'head=',h.temp[which(h.temp<p.min.surf)]))
			write(file=logfile,append=T,x=paste('SIMULATION STOPPED',date()))
			update.data()
			stop_() 
		}
		if (any(h.temp>p.max.surf)) {
                        write(file=logfile,append=T,paste('high head in ponds:',p.names[which(h.temp>p.max.surf)]))
                        write(file=logfile,append=T,paste('max head=',
				p.max.surf[which(h.temp>p.max.surf)],'head=',h.temp[which(h.temp>p.max.surf)]))
			write(file=logfile,append=T,x=paste('SIMULATION STOPPED',date()))
			update.data()
                        stop_()
                }		

		#Check for high flows#######################################################
		if (sum(ff[8:16,t])>=3*mgd.to.cfs) {
			write(file=logfile,append=T,x=paste('Flow to marsh box > 3 mgd, must turn on pond pumps',date()))
			write(file=logfile,append=T,x=paste('SIMULATION STOPPED',date()))
			update.data()
			stop_()
		}
	}
    write(file=paste(path,'modelRunning',sep=''),'1')
